####################################################################################
 # Replication file for "Local Ownership of IMF Conditionality Programs"
 #
 # Nikitas Konstantinidis and Bernhard Reinsberg
 # Contact: bernhard.reinsberg@glasgow.ac.uk
####################################################################################

# Estimating an ownership measure 
####################################################################################

# Data 

  library(foreign)
  library(DataCombine)
  library(Synth)
  library(stargazer)
  library(dplyr)

  setwd("C:\\Berni\\Forschung\\IMFOWN\\Submissions\\ISQ")

  #d=read.dta("IMF spells.dta")
  d=read.dta("Unimputed data.dta")
  #D=read.dta("IMF data.dta")
  D=read.dta("Imputed data.dta")
  own=read.dta("TSCS_by_sector.dta")
  own=own[which(own$failed==0 | own$failed==2), ]


# Definition of treatment sets and control sets (obeying cutoff rules)  
  
  Spells=which(d$spell==1)
  Treated.EXT=Control.FIN=Control.Any=numeric(0)
  Control.EXT=Treated.FIN=Treated.Any=Treated=Spells[1]
  
  ## 5-year gaps between spells ##  
  MinGap=5
  for (s in 2:length(Spells)){
    diff=ifelse(d$cid[Spells[s]]==d$cid[Spells[s-1]], d$year[Spells[s]]-d$year[Spells[s-1]], 99)
    if (diff>=MinGap){
      Treated=c(Treated, Spells[s])
    }
    
    p.obs=ifelse(s<length(Spells), 
                 (Spells[s]:(Spells[s+1]-1))[d$IMFnn[Spells[s]:(Spells[s+1]-1)]==1],
                 Spells[s])
    
    if (any(d$SCsEXT[p.obs], na.rm=T)>0){
      Treated.EXT=c(Treated.EXT, Spells[s])
    } else {
      Control.EXT=c(Control.EXT, Spells[s])
    }
    
    if (any(d$SCsFIN[p.obs], na.rm=T)>0){
      Treated.FIN=c(Treated.FIN, Spells[s])
    } else {
      Control.FIN=c(Control.FIN, Spells[s])
    }
    
    if (any(d$SCsEXT[p.obs]+d$SCsFIN[p.obs], na.rm=T)>0){
      Treated.Any=c(Treated.Any, Spells[s])
    } else {
      Control.Any=c(Control.Any, Spells[s])
    }
  }
  
  Treated.EXT=intersect(Treated, Treated.EXT)
  Control.EXT=intersect(Treated, Control.EXT)
  Treated.FIN=intersect(Treated, Treated.FIN)
  Control.FIN=intersect(Treated, Control.FIN)  
  Treated.Any=intersect(Treated, Treated.Any)
  Control.Any=intersect(Treated, Control.Any)
  
  length(Treated.EXT)
  length(Treated.FIN)
  length(Treated.Any)
  
  ## Statistics on sample selection ##
  length(Spells)
  length(Treated)
  length(Treated.EXT)
  length(Treated.FIN)
  library(cobalt)
  D$treat=0
  D$treat[Spells]=1
  covs=subset(D, select=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp","BA2TOT","BA2scope","loanquota","nUnder","treat"))
  bal.tab(covs, treat=D$treat, disp=c("means","sds"))
  bal.tab.t=data.frame(covs=names(covs), p=NA)
  for (j in 1:(ncol(covs)-1)){
    bal.tab.t$p[j]=t.test(covs[covs$treat==1,j], covs[covs$treat==0,j])$p.value
  }
  (bal.tab.t)
  D$treat=NA
  D$treat[Spells]=0
  D$treat[Treated]=1
  covs=subset(D[is.na(D$treat)==F,], select=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp","BA2TOT","BA2scope","loanquota","nUnder","treat"))
  bal.tab(covs, treat=D$treat[is.na(D$treat)==F], disp=c("means","sds"))
  bal.tab.t=data.frame(covs=names(covs), p=NA)
  for (j in 1:(ncol(covs)-1)){
    bal.tab.t$p[j]=t.test(covs[covs$treat==1,j], covs[covs$treat==0,j])$p.value
  }
  (bal.tab.t)
  D$treat=NA
  D$treat[Treated]=0
  D$treat[Treated.Any]=1
  covs=subset(D[is.na(D$treat)==F,], select=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp","BA2TOT","BA2scope","loanquota","nUnder","treat"))
  bal.tab(covs, treat=D$treat[is.na(D$treat)==F], disp=c("means","sds"))
  bal.tab.t=data.frame(covs=names(covs), p=NA)
  for (j in 1:(ncol(covs)-1)){
    bal.tab.t$p[j]=t.test(covs[covs$treat==1,j], covs[covs$treat==0,j])$p.value
  }
  (bal.tab.t)
  D$treat=NA
  D$treat[Treated.Any]=0
  U=unique(own$spellid)
  D$treat[Treated.Any[U]]=1
  covs=subset(D[is.na(D$treat)==F,], select=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp","BA2TOT","BA2scope","loanquota","nUnder","treat"))
  bal.tab(covs, treat=D$treat[is.na(D$treat)==F], disp=c("means","sds"))
  bal.tab.t=data.frame(covs=names(covs), p=NA)
  for (j in 1:(ncol(covs)-1)){
    bal.tab.t$p[j]=t.test(covs[covs$treat==1,j], covs[covs$treat==0,j])$p.value
  }
  (bal.tab.t)
  
  
  # Propensity score for conditionality treatment ##
  PS.EXT=glm(treatedEXT~nUnder+BA2TOT+BA2scope+ s_unga3g7+unsc_DSV+lngdppc+gdp_growth_WDI+cab_gdp_WDI+debtser_gni_WDI+reserves_WDI+p_polity2_QOG+as.factor(regid)+as.factor(incid)+as.factor(year), data=D, family="binomial")
  D$PS.EXT=predict(PS.EXT, type="response")
  PS.FIN=glm(treatedFIN~nUnder+BA2TOT+BA2scope+ s_unga3g7+unsc_DSV+lngdppc+gdp_growth_WDI+cab_gdp_WDI+debtser_gni_WDI+reserves_WDI+p_polity2_QOG+as.factor(regid)+as.factor(incid)+as.factor(year), data=D, family="binomial")
  D$PS.FIN=predict(PS.FIN, type="response")
  PS.EXT=PS.FIN=NULL
  
  set.seed(13)
  
  iterations=7
  se=1

# User-written functions

  ## Function to shift start years ##
  CreateControlData=function(cid, year, sector){
    obs.time=(max(1980,year-6)):(min(2014,year+6))
    dat=D[D$cid==cid & D$year%in%obs.time,]
    dat$cid=0
    dat$Name=D$iso3[D$cid==cid & D$year==year]
  
    Control=intersect(Control.EXT, Control.FIN)
  
    for (k in 1:length(Control)){
      k.cid=d$cid[Control[k]]
      k.year=d$year[Control[k]]  
    
    if (k.year>1980 & k.year<2014){
      new.dat=D[D$cid==k.cid & D$year>=max(1980,k.year-6) & D$year<=min(2014,k.year+6),]
      new.dat$year[new.dat$spell==1]=year
      after=(which(new.dat$spell==1)+1):nrow(new.dat)
      new.dat$year[after]=seq(year+1,length.out=length(after),by=1)
      before=1:(which(new.dat$spell==1)-1)
      new.dat$year[before]=sort(seq(year-1,length.out=length(before),by=-1))
      new.dat$cid=k
      iso3=new.dat$iso3[1]
      merge.dat=expand.grid(cid=k, year=obs.time)
      merge.dat=merge(merge.dat, new.dat, by=c("cid","year"), all.x=T, all.y=F)
      merge.dat$Name=paste(iso3,k,sep="")
      rows.missing=which(is.na(merge.dat[,4]))
      rows.exist=which(is.na(merge.dat[,4])==F)
      if (min(rows.exist,na.rm=T)>max(rows.missing,na.rm=T)){
        merge.dat[rows.missing,3:ncol(merge.dat)]=merge.dat[min(rows.exist),3:ncol(merge.dat)]
      } else {
        merge.dat[rows.missing,3:ncol(merge.dat)]=merge.dat[max(rows.exist),3:ncol(merge.dat)]
      }
      dat=rbind(dat, merge.dat)
    }
  }
  return=dat
  }

  ## Function to optimize plotting window ##
  Lim=function(v1,v2,dir){
    if (dir==0){
      v=min(c(v1,v2),na.rm=T)
      v=ifelse(v>0, 0.64*v, 1.36*v)
    }
    else {
      v=max(c(v1,v2),na.rm=T)
      v=ifelse(v>0, 1.36*v, 0.64*v)
    }
    return=v
  }

# Synthetic control method

  ## EXT sector ##

  synthSE.TR=matrix(NA, ncol=length(Treated.EXT), nrow=13)
  synthSE.CR=synthSE.TR
  synthSE.se=synthSE.TR
  synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
  synthSE.diff=expand.grid(iso3=d$iso3[Treated.EXT], t=0:6, diff=NA)
  synthSE.weights.var=matrix(NA, ncol=length(Treated.EXT), nrow=18+2)

  pdf("EXT-EcGI-SE.pdf")

for (i in 1:length(Treated.EXT)){
  
  treated=Treated.EXT 
  control=Control.EXT 
  Y1="KOFEcGIdj"
  Y0="trade_WDI"
  X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  cty=d$cid[treated[i]]
  yr=d$year[treated[i]]
  D.i=CreateControlData(cty,yr,"Any")
  
  ## Matching
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  
  ## Main routine  
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  
    y0.out=data.frame(synth.tables$tab.w[synth.tables$tab.w$w.weights>0,])
    y0.out=cbind(iso3=D.i$iso3[1], year=d$year[Treated.EXT[i]], y0.out)
    write.table(y0.out, "Ownership-cases-EXT.txt", append=T)
    
    v.out=data.frame(vars=rownames(synth.tables$tab.v), weight=as.numeric(synth.tables$tab.v)) 
    v.out=cbind(iso3=D.i$iso3[1], year=d$year[Treated.EXT[i]], v.out)
    write.table(v.out, "Ownership-variables-EXT.txt", append=T)
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Main=paste(d$cname[treated[i]],yr), Ylab="De-jure economic liberalization", Xlab=paste("|D|=",length(ctrls.cid),sep=""), 
            Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[i,1]<-d$iso3[treated[i]]
  synthSE.diag[i,2]<-yr
  synthSE.diag[i,3]<-pre.years 
  synthSE.diag[i,4]<-post.years 
  synthSE.diag[i,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[i,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[i,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  synthSE.diff$diff[(synthSE.diff$iso3 %in% d$iso3[treated[i]])[1:(post.years+1)]]<-post.diff
  synthSE.weights.var[1:18,i]=as.numeric(synth.tables$tab.v)
  synthSE.weights.var[19,i]=mahalanobis(synth.tables$tab.pred[1:12,1],synth.tables$tab.pred[1:12,2],cov=cov(D[,X1]))
    M=data.frame(T=synth.tables$tab.pred[1:12,1], SC=synth.tables$tab.pred[1:12,2], sd=sqrt(diag(cov(D[,X1]))), w=as.numeric(synth.tables$tab.v)[1:12])
    M$t=(M$T-M$SC)/M$sd
  synthSE.weights.var[20,i]=weighted.mean(M$t, M$w)
  
  ## Adding uncertainty estimates 
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,i]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,i]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,i]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,i]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,i]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,i]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,i]<-dat.i.k$Y1plot
      synthSE.CR[,i]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,i]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
  }
}

  dev.off()
  write.csv(synthSE.diag, "EXT-EcGI-SE_1.csv")
  write.csv(synthSE.TR, "EXT-EcGI-SE_2.csv")
  write.csv(synthSE.CR, "EXT-EcGI-SE_3.csv")
  write.csv(synthSE.se, "EXT-EcGI-SE_4.csv")  
  write.csv(synthSE.diff, "EXT-EcGI-SE_5.csv")
  write.csv(synthSE.weights.var, "EXT-EcGI-w-EXT.csv")

  save.image(file="R-EXT-EcGI.Rdata")

  
  ## FIN sector ##

  synthSE.TR=matrix(NA, ncol=length(Treated.FIN), nrow=13)
  synthSE.CR=synthSE.TR
  synthSE.se=synthSE.TR
  synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
  synthSE.diff=expand.grid(iso3=d$iso3[Treated.FIN], t=0:6, diff=NA)
  synthSE.weights.var=matrix(NA, ncol=length(Treated.FIN), nrow=18+2)

  pdf("FIN-FiGI-SE.pdf")

for (i in 1:length(Treated.FIN)){
  
  treated=Treated.FIN
  control=Control.FIN 
  Y1="KOFFiGIdj"
  Y0="fdi_gdp_WDI"
  X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  cty=d$cid[treated[i]]
  yr=d$year[treated[i]]
  D.i=CreateControlData(cty,yr,"Any")
  
  ## Matching
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  
  ## Main routine  
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  
    y0.out=data.frame(synth.tables$tab.w[synth.tables$tab.w$w.weights>0,])
    y0.out=cbind(iso3=D.i$iso3[1], year=d$year[Treated.FIN[i]], y0.out)
    write.table(y0.out, "Ownership-cases-FIN.txt", append=T)
  
    v.out=data.frame(vars=rownames(synth.tables$tab.v), weight=as.numeric(synth.tables$tab.v)) 
    v.out=cbind(iso3=D.i$iso3[1], year=d$year[Treated.FIN[i]], v.out)
    write.table(v.out, "Ownership-variables-FIN.txt", append=T)
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Main=paste(d$cname[treated[i]],yr), Ylab="De-jure financial liberalization", Xlab=paste("|D|=",length(ctrls.cid),sep=""), 
            Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[i,1]<-d$iso3[treated[i]]
  synthSE.diag[i,2]<-yr
  synthSE.diag[i,3]<-pre.years 
  synthSE.diag[i,4]<-post.years 
  synthSE.diag[i,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[i,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[i,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  synthSE.diff$diff[(synthSE.diff$iso3 %in% d$iso3[treated[i]])[1:(post.years+1)]]<-post.diff
  synthSE.weights.var[1:18,i]=as.numeric(synth.tables$tab.v)
  synthSE.weights.var[19,i]=mahalanobis(synth.tables$tab.pred[1:12,1],synth.tables$tab.pred[1:12,2],cov=cov(D[,X1]))
    M=data.frame(T=synth.tables$tab.pred[1:12,1], SC=synth.tables$tab.pred[1:12,2], sd=sqrt(diag(cov(D[,X1]))), w=as.numeric(synth.tables$tab.v)[1:12])
    M$t=(M$T-M$SC)/M$sd
  synthSE.weights.var[20,i]=weighted.mean(M$t, M$w)
  
  ## Adding uncertainty estimates 
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,i]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,i]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,i]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,i]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,i]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,i]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,i]<-dat.i.k$Y1plot
      synthSE.CR[,i]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,i]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
  }
}

  dev.off()
  write.csv(synthSE.diag, "FIN-FiGI-SE_1.csv")
  write.csv(synthSE.TR, "FIN-FiGI-SE_2.csv")
  write.csv(synthSE.CR, "FIN-FiGI-SE_3.csv")
  write.csv(synthSE.se, "FIN-FiGI-SE_4.csv")  
  write.csv(synthSE.diff, "FIN-FiGI-SE_5.csv")
  write.csv(synthSE.weights.var, "FIN-FiGI-w.csv")

  save.image(file="R-FIN-FiGI.Rdata")

  
  
# Diagnostics with the ownership measure
####################################################################################

# Load SCM output with annual ownership and goodness-of-fit ownership  
  
  #d=read.dta("IMF spells.dta")
  d=read.dta("Unimputed data.dta")
  
    ## Are sectoral conditions representative of wider demands?
    d$SCsboth=d$SCsEXT+d$SCsFIN
    imf=which(d$IMFnn==1)
    cor(d$SCsboth[imf], d$SCB[imf], use="pairwise.complete.obs")
    cor(d$SCsEXT[imf], d$SCB[imf], use="pairwise.complete.obs")
    cor(d$SCsFIN[imf], d$SCB[imf], use="pairwise.complete.obs")
    pairs(d[imf,c("SCsboth","SCB","QCB","BA2TOT","BA2scope")])
    pairs(d[imf,c("SCsEXT","SCB","QCB","BA2TOT","BA2scope")])
    pairs(d[imf,c("SCsFIN","SCB","QCB","BA2TOT","BA2scope")])
    --> low correlation (if not negative)
  
  own=read.dta("TSCS_by_sector.dta")
  own=own[which(own$failed==0 | own$failed==2), ]
  ## exclude spells with any interruption in t1, t2, t3, or t4 (==0)
  ## still include negative deviations where no implementation information was available (==2)
  
  rmspe=read.dta("RMSPE_by_sector.dta")
  rmspe=rmspe[which(rmspe$errorcode==0 | rmspe$errorcode==2), ]
  
  rmspe=merge(rmspe, d, by=c("iso3","year"), all.x=T, all.y=F, suffixes=c("","2"))
  rmspe$SCsboth=rmspe$SCsEXT+rmspe$SCsFIN
  own=merge(own, d, by=c("iso3","year"), all.x=T, all.y=F, suffixes=c("","2"))
  own$SCsboth=own$SCsEXT+own$SCsFIN
  
  rmspe$fut5=rmspe$w=NA
  for (i in 1:nrow(rmspe)){
    cty=rmspe$iso3[i]
    year=rmspe$year[i]
    rmspe$fut5[i]=mean(d$IMFnn[which(d$iso3%in%cty & d$year%in%(year+4):(year+8))], na.rm=T)
    rmspe$w[i]=rmspe$BA2TOT[i]-rmspe$cBATOT[i]
  }
  
  own$w=own$sw=NA
  own$fut5=NA
  for (i in 1:nrow(own)){
    idx=which(d$iso3%in%own$iso3[i] & d$year%in%own$year[i])
    if (length(idx)>0){
      own$w[i]=(d$BA2TOT[idx]-d$cBATOT[idx])
      own$sw[i]=own$w[i]/d$BA2TOT[idx]
    }  
    after=(own$year[i]+4):(own$year[i]+8)
    own$fut5[i]=mean(d$IMFnn[which(d$iso3%in%own$iso3[i] & d$year%in%after)], na.rm=T)
  }
  
  #NOT - embargoed data
  #revtscs=read.dta("E:\\_imf\\GHFPRT\\Data\\IMF_reviews_data.dta")
  #own$inttemp=own$intpermcor=own$adur=NULL
  #own=merge(own, revtscs, by=c("iso3","year"), all.x=T, all.y=F)
  
  
  ## Descriptive information on ownership 
  rmspe$streamlining=rmspe$year>2000
  data.frame(rmspe %>% group_by(streamlining) %>% summarise(y0=mean(o), y1=sd(o), y2=min(o), y3=max(o)))
  data.frame(rmspe %>% group_by(region_WDI) %>% summarise(y0=mean(o), y1=sd(o), y2=min(o), y3=max(o)))
  
  ## Descriptive statistics for supplemental appendix 
  stargazer(d[,c("SCsEXT","SCsFIN","KOFEcGIdj","KOFEcGIdf")], type="text", summary=T)
  
  
# Diagnostic plots 
  
  ## Determinants
  pdf("A1-o-RMSPE.pdf")
  plot(rmspe$SCsboth[], rmspe$oT[], xlab="Structural conditions", ylab="Ownership")
  lm0=lm(oT~SCsboth, data=rmspe[,])
  lines(lm0$model$SCsboth, lm0$fitted.values, lty=2, col="gray")
  rug(jitter(rmspe$SCsboth, 0.08), col="blue")
  ## Figure 2 ## 
  dev.off() 
  
    ## remove outlier
    pdf("Figure R1.pdf")
    L=which(rmspe$SCsboth>20)
    plot(rmspe$SCsboth[-L], rmspe$oT[-L], xlab="Structural conditions", ylab="Ownership")
    lm0=lm(oT~SCsboth, data=rmspe[-L,])
    lines(lm0$model$SCsboth, lm0$fitted.values, lty=2, col="gray")
    rug(jitter(rmspe$SCsboth, 0.08), col="blue")
    dev.off()
    
    ## LOO-CV test
    X=data.frame(coef=numeric(nrow(rmspe)), t=rep(NA, nrow(rmspe)))
    for (j in 1:nrow(rmspe)){
      CV=setdiff(1:nrow(rmspe), j)
      lm.CV=summary(lm(oT~SCsboth, data=rmspe[CV,]))
      X[j,1]=lm.CV$coefficients[2,1]
      X[j,2]=abs(lm.CV$coefficients[2,3])
    }
    apply(X,2,mean)
    summary(X)
    sqrt(diag(var(X)))
    
  
  #S.EXT=which(rmspe$sectorid==1)
  #pdf("A1-o-EXT.pdf")
  #plot(rmspe$SCsEXT[S.EXT], rmspe$oT[S.EXT], xlab="External sector conditions", ylab="Ownership")
  #lm0=lm(oT~SCsEXT, data=rmspe[S.EXT,])
  #lines(lm0$model$SCsEXT, lm0$fitted.values, lty=2, col="gray")
  #rug(jitter(rmspe$SCsEXT[S.EXT], 0.08), col="blue")
  #dev.off() 
  
  #S.FIN=which(rmspe$sectorid==2)
  #pdf("A1-o-FIN.pdf")
  #plot(rmspe$SCsFIN[S.FIN], rmspe$oT[S.FIN], xlab="Financial sector conditions", ylab="Ownership")
  #lm0=lm(oT~SCsFIN, data=rmspe[S.FIN,])
  #lines(lm0$model$SCsFIN, lm0$fitted.values, lty=2, col="gray")
  #rug(jitter(rmspe$SCsFIN[S.FIN], 0.08), col="blue")
  #dev.off() 
  
  pdf("A2-o-RMSPE.pdf")
  plot(rmspe$PAsTOT, rmspe$oT[], xlab="Prior actions", ylab="Ownership")
  lm1=lm(oT~PAsTOT, data=rmspe[,])
  lines(lm1$model$PAsTOT, lm1$fitted.values, lty=2, col="gray")
  rug(jitter(rmspe$PAsTOT[], 0.08), col="blue")
  ## Figure 3 ## 
  dev.off() 
  
    ## remove outlier
    pdf("Figure R2.pdf")
    L=which(rmspe$PAsTOT>12)
    plot(rmspe$PAsTOT[-L], rmspe$oT[-L], xlab="Prior actions", ylab="Ownership")
    lm1=lm(oT~PAsTOT, data=rmspe[-L,])
    lines(lm1$model$PAsTOT, lm1$fitted.values, lty=2, col="gray")
    rug(jitter(rmspe$PAsTOT[-L], 0.08), col="blue")
    dev.off()
    
    ## LOO-CV test
    X=data.frame(coef=numeric(nrow(rmspe)), t=rep(NA, nrow(rmspe)))
    for (j in 1:nrow(rmspe)){
      CV=setdiff(1:nrow(rmspe), j)
      lm.CV=summary(lm(oT~PAsTOT, data=rmspe[CV,]))
      X[j,1]=lm.CV$coefficients[2,1]
      X[j,2]=abs(lm.CV$coefficients[2,3])
    }
    apply(X,2,mean)
    summary(X)
    sqrt(diag(var(X)))
  
  
  #pdf("A2-o-EXT.pdf")
  #plot(rmspe$PAsTOT[rmspe$sectorid==1], rmspe$oT[rmspe$sectorid==1], xlab="Prior actions", ylab="Ownership")
  #lm1=lm(oT~PAsTOT, data=rmspe[rmspe$sectorid==1,])
  #lines(lm1$model$PAsTOT, lm1$fitted.values, lty=2, col="gray")
  #rug(jitter(rmspe$PAsTOT[rmspe$sectorid==1], 0.08), col="blue")
  #dev.off() 
  
  #pdf("A2-o-FIN.pdf")
  #plot(rmspe$PAsTOT[S.FIN], rmspe$oT[S.FIN], xlab="Prior actions", ylab="Ownership")
  #lm1=lm(oT~PAsTOT, data=rmspe[S.FIN,])
  #lines(lm1$model$PAsTOT, lm1$fitted.values, lty=2, col="gray")
  #rug(jitter(rmspe$PAsTOT[S.FIN], 0.08), col="blue")
  #dev.off() 
  
  
  pdf("B2-o-RMSPE.pdf")
  S=which(rmspe$fut5<1)
  plot(rmspe$oT[S], rmspe$fut5[S], xlab="Ownership", ylab="Post-program participation rate", main=" ")
  lmi=lm(fut5~oT, data=rmspe[S,])
  lines(lmi$model$oT, lmi$fitted.values, lty=2, col="gray") 
  rug(jitter(rmspe$oT[S], 0.08), col="blue")
  ## Figure 4 ## 
  dev.off()
  
  #pdf("B2-o-EXT.pdf")
  #S.o.EXT=which(rmspe$sectorid==1 & rmspe$fut5<1)
  #plot(rmspe$oT[S.o.EXT], rmspe$fut5[S.o.EXT], xlab="Ownership", ylab="Post-program participation rate", main=" ")
  #lmi=lm(fut5~oT, data=rmspe[S.o.EXT,])
  #lines(lmi$model$oT, lmi$fitted.values, lty=2, col="gray") 
  #rug(jitter(rmspe$oT[S.o.EXT], 0.08), col="blue")
  #dev.off()
  
  #pdf("B2-o-FIN.pdf")
  #S.o.FIN=which(rmspe$sectorid==2 & rmspe$fut5<1)
  #plot(rmspe$oT[S.o.FIN], rmspe$fut5[S.o.FIN], xlab="Ownership", ylab="Post-program participation rate", main=" ")
  #lmi=lm(fut5~oT, data=rmspe[S.o.FIN,])
  #lines(lmi$model$oT, lmi$fitted.values, lty=2, col="gray") 
  #rug(jitter(rmspe$oT[S.o.FIN], 0.08), col="blue")
  #dev.off()
  
  
  rmspe$dj1=rmspe$dj2=rmspe$df1=rmspe$df2=rmspe$trdfmdj=rmspe$fidfmdj=rmspe$DD=rmspe$DDeco=NA
  for (i in 1:nrow(rmspe)){
    cty=rmspe$iso3[i]
    year=rmspe$year[i]
    idx1=which(d$iso3%in%cty & d$year%in%year)
    idx2=which(d$iso3%in%cty & d$year%in%(year+3))
    
    if (rmspe$sectorid[i]==1){
      rmspe$dj1[i]=d$KOFTrGIdj[idx1]
      rmspe$dj2[i]=ifelse(length(idx2)>0, d$KOFTrGIdj[idx2], NA)
      rmspe$df1[i]=d$KOFTrGIdf[idx1]
      rmspe$df2[i]=ifelse(length(idx2)>0, d$KOFTrGIdf[idx2], NA)
      rmspe$trdfmdj[i]=(rmspe$df2[i]-rmspe$df1[i])-(rmspe$dj2[i]-rmspe$dj1[i])
    }
    
    if (rmspe$sectorid[i]==2){
      rmspe$dj1[i]=d$KOFFiGIdj[idx1]
      rmspe$dj2[i]=ifelse(length(idx2)>0, d$KOFFiGIdj[idx2], NA)
      rmspe$df1[i]=d$KOFFiGIdf[idx1]
      rmspe$df2[i]=ifelse(length(idx2)>0, d$KOFFiGIdf[idx2], NA)
      rmspe$fidfmdj[i]=(rmspe$df2[i]-rmspe$df1[i])-(rmspe$dj2[i]-rmspe$dj1[i])
    }    
    
    rmspe$DD[i]=ifelse(rmspe$sectorid[i]==1, rmspe$trdfmdj[i], rmspe$fidfmdj[i])
    rmspe$DDeco[i]=ifelse(length(idx2)>0, (d$KOFEcGIdf[idx2]-d$KOFEcGIdf[idx1])-(d$KOFEcGIdj[idx2]-d$KOFEcGIdj[idx1]), NA)
  }  
  
  ## Double difference
  pdf("C1-o-RMSPE.pdf")
  plot(rmspe$oT[], rmspe$DD[], xlab="Ownership", ylab="De-facto -- de-jure")
  lmdd=lm(DD~oT, data=rmspe[,])
  lines(lmdd$model$oT, lmdd$fitted.values, lty=2, col="gray")
  rug(jitter(rmspe$oT[], 0.08), col="blue")
  ## Figure 5 ##
  dev.off()
  
  #pdf("C1-o-EXT.pdf")
  #plot(rmspe$oT[S.EXT], rmspe$trdfmdj[S.EXT], xlab="Ownership", ylab="De-facto -- de-jure")
  #lmdd=lm(trdfmdj~oT, data=rmspe[S.EXT,])
  #lines(lmdd$model$oT, lmdd$fitted.values, lty=2, col="gray")
  #rug(jitter(rmspe$rmspe[S.EXT], 0.08), col="blue")
  #dev.off()
  
  #pdf("C1-o-FIN.pdf")
  #plot(rmspe$oT[S.FIN], rmspe$fidfmdj[S.FIN], xlab="Ownership", ylab="De-facto -- de-jure")
  #lmdd=lm(fidfmdj~oT, data=rmspe[S.o.FIN,])
  #lines(lmdd$model$oT, lmdd$fitted.values, lty=2, col="gray")
  #rug(jitter(rmspe$oT[S.FIN], 0.08), col="blue")
  #dev.off()
  
  rmspe$waiv2=c(4, 5, 14, 13, 13, 5, NA, 0, 0, 12, NA, NA, NA, NA, 5, 8, NA, 11, 11, NA, 1, 1, NA, 1, NA, 20, 20, 6, 6, 0, 0, 0)
  ## Data on IMF implementation shared in confidence: waiver over entire program length ## 
  pdf("B3-o-RMSPE.pdf")
  plot(rmspe$oT[], rmspe$waiv2[], xlab="Ownership", ylab="Number of waived conditions", main=" ")
  lmi=lm(waiv2~oT, data=rmspe[,])
  lines(lmi$model$oT, lmi$fitted.values, lty=2, col="gray") 
  rug(jitter(rmspe$rmspe[S], 0.08), col="blue")
  ## Figure A7 - unused ##
  dev.off()
  
  pdf("B4-o-RMSPE.pdf")
  plot(rmspe$oT[], rmspe$w[], xlab="Ownership", ylab="Waived conditions at entry", main=" ")
  lmi=lm(w~oT, data=rmspe[,])
  lines(lmi$model$oT, lmi$fitted.values, lty=2, col="gray") 
  rug(jitter(rmspe$rmspe[], 0.08), col="blue")
  ## Figure A7 ##
  dev.off()
  
  
  #pdf("B3-o-EXT.pdf")
  #plot(rmspe$oT[S.EXT], rmspe$waiv2[S.EXT], xlab="Ownership", ylab="Number of waived conditions", main=" ")
  #lmi=lm(waiv2~oT, data=rmspe[S.EXT,])
  #lines(lmi$model$oT, lmi$fitted.values, lty=2, col="gray") 
  #rug(jitter(rmspe$rmspe[S.EXT], 0.08), col="blue")
  #dev.off()
  
  #pdf("B3-o-FIN.pdf")
  #plot(rmspe$oT[S.FIN], rmspe$waiv2[S.FIN], xlab="Ownership", ylab="Number of waived conditions", main=" ")
  #lmi=lm(waiv2~oT, data=rmspe[S.FIN,])
  #lines(lmi$model$oT, lmi$fitted.values, lty=2, col="gray") 
  #rug(jitter(rmspe$rmspe[S.FIN], 0.08), col="blue")
  #dev.off()
  
  pdf("A3-o-RMSPE.pdf")
  plot(rmspe$fi_index_cl_QOG, rmspe$oT, xlab="Economic freedom index", ylab="Ownership", ylim=c(-0.12,0))  
  lm7=lm(oT~fi_index_cl_QOG, data=rmspe[,])
  lines(lm7$model$fi_index_cl_QOG, lm7$fitted.values, lty=2, col="gray")
  rug(jitter(own$fi_index_cl_QOG, 0.08), col="blue")
  ## Figure A4 ## 
  dev.off()
  
  pdf("A4-o-RMSPE.pdf")
  plot(rmspe$foreignbkct_cv[], rmspe$oT[], xlab="Number of foreign banks", ylab="Ownership")  
  lm11=lm(oT~foreignbkct_cv, data=rmspe[,])
  lines(lm11$model$foreignbkct_cv, lm11$fitted.values, lty=2, col="gray")
  rug(jitter(rmspe$foreignbkct_cv[], 0.08), col="blue")
  ## Figure A5 ##
  dev.off()
  
  rmspe$lnbkus=log(1+rmspe$bkexp_repUSA)  
  pdf("A4-o-US.pdf")
  plot(rmspe$lnbkus[], rmspe$oT[], xlab="Logged exposure by US banks", ylab="Ownership")  
  lm12=lm(oT~lnbkus, data=rmspe[,])
  lines(lm12$model$lnbkus, lm12$fitted.values, lty=2, col="gray")
  rug(jitter(rmspe$lnbkus[], 0.08), col="blue")
  ## Figure A5 - unused ##
  dev.off() 
  
  
  
# Supplemental appendix: Additional validation plots
####################################################
  
  # Annual
  #T=which(own$t==0)
  
  #pdf("A1-o-TSCS.pdf")
  #plot(own$SCsboth[T], own$oT3[T], xlab="Structural conditions", ylab="Ownership", sub="t=0")
  #lm0=lm(oT3~SCsboth, data=own[T,])
  #lines(lm0$model$SCsboth, lm0$fitted.values, lty=2, col="gray")
  #rug(jitter(own$SCsboth, 0.08), col="blue")
  #dev.off() 
  
  y0=which(own$t==0)
  y1=which(own$t==1)
  y2=which(own$t==2)
  y3=which(own$t==3)
  
  pdf("A1-o-all.pdf")
  plot(own$SCsboth[y0], own$oT3[y0], xlab="Structural conditions", ylab="Ownership")
  points(own$SCsboth[y1], own$oT3[y1], col=3)
  points(own$SCsboth[y2], own$oT3[y2], col=5)
  points(own$SCsboth[y3], own$oT3[y3], col=6)
  lm1=lm(oT3~SCsboth, data=own[y0,])
  lines(lm1$model$SCsboth, lm1$fitted.values, lty=4, col="gray")
  legend("bottomright", bg="white", legend=c("t=0","t=1","t=2","t=3"),col=c(1,3,5,6), pch=c(1,1,1,1), cex=0.8)  
  ## Figure A2 ##
  dev.off()
  
  
  pdf("A2-o-TSCS.pdf")
  plot(own$PAsTOT[y0], own$oT3[y0], xlab="Prior actions", ylab="Ownership")
  lm1=lm(oT3~PAsTOT, data=own[y0,])
  lines(lm1$model$PAsTOT, lm1$fitted.values, lty=2, col="gray")
  rug(jitter(own$PAsTOT[y0], 0.08), col="blue")
  ## Figure A3 ##
  dev.off() 
  
  
  pdf("B2-o-TSCS.pdf")
  plot(own$oT3[], own$fut5[], xlab="Ownership", ylab="Post-program participation rate", main=" ")
  lmi=lm(fut5~oT3, data=own[,])
  lines(lmi$model$oT3, lmi$fitted.values, lty=2, col="gray") 
  rug(jitter(own$diff[], 0.08), col="blue")
  ## Figure A6 ##
  dev.off()
  
  pdf("B3-o-TSCS.pdf")
  plot(own$oT3[], own$w[], xlab="Ownership", ylab="Number of waived conditions", main=" ")
  lmi=lm(w~oT3, data=own[,])
  lines(lmi$model$oT3, lmi$fitted.values, lty=2, col="gray") 
  rug(jitter(own$diff[], 0.08), col="blue")
  ## Figure A7 - not shown ##
  dev.off()
  
  
# Case-based robustness of SCM approach using 'IDN 1997'
########################################################
  
  se=1
  iterations=7
  
  synthSE.TR=matrix(NA, ncol=1, nrow=13)
  synthSE.CR=synthSE.TR
  synthSE.se=synthSE.TR
  synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
  synthSE.diff=expand.grid(iso3="IDN", t=0:6, diff=NA)
  
  d$iso3[Treated.EXT]
  idn=16
  i=idn
  treated=Treated.EXT 
  control=Control.EXT 
  Y1="KOFEcGIdj"
  Y0="trade_WDI"
  X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  
# Reproduce output from the main routine (for 'IDN 1997')  

    cty=d$cid[treated[i]]
    yr=d$year[treated[i]]
    D.i=CreateControlData(cty,yr,"Any")
    
    ## Matching
    
    ctrls.cid=setdiff(unique(D.i$cid), 0)
    
    ## Main routine  
    
    dat.i=dataprep(D.i,
                   predictors=X1, 
                   predictors.op="mean",
                   dependent=Y1,
                   unit.variable="cid",
                   time.variable="year",
                   special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                           list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                           list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                   treatment.identifier=0,
                   controls.identifier=ctrls.cid,
                   time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                   time.optimize.ssr=(max(1980,yr-6):yr),
                   unit.names.variable="Name",
                   time.plot=(max(1980,yr-6):min(2014,yr+6))
    )
    
    if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
      idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
      dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
    }
    
    synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
    if (inherits(synth.out, "error")) next
    
    synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
    print(synth.tables)
    
    y0.out=data.frame(synth.tables$tab.w[synth.tables$tab.w$w.weights>0,])
    y0.out=cbind(iso3=D.i$iso3[1], year=d$year[Treated.EXT[i]], y0.out)
    y0.out
    
    v.out=data.frame(vars=rownames(synth.tables$tab.v), weight=as.numeric(synth.tables$tab.v)) 
    v.out=cbind(iso3=D.i$iso3[1], year=d$year[Treated.EXT[i]], v.out)
    v.out
    
    M=data.frame(T=synth.tables$tab.pred[1:12,1], SC=synth.tables$tab.pred[1:12,2], sd=sqrt(diag(cov(D[,X1]))), w=as.numeric(synth.tables$tab.v)[1:12])
    M$t=(M$T-M$SC)/M$sd
    (M)
    (weighted.mean(M$t, M$w))
    
    Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
    
    pdf("IDN-1997.pdf")
    path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
    legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
    abline(v=yr, lty=2, col="gray")
    abline(h=0, lty=2, col="gray")
    
    pre.years=6+(min(yr,1986)-1986)
    post.years=6-(max(yr,2008)-2008)
    post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
    pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
    
    synthSE.diag[i,1]<-d$iso3[treated[i]]
    synthSE.diag[i,2]<-yr
    synthSE.diag[i,3]<-pre.years 
    synthSE.diag[i,4]<-post.years 
    synthSE.diag[i,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
    synthSE.diag[i,6]<-sqrt(mean(pre.diff^2,na.rm=T))
    synthSE.diag[i,7]<-sqrt(mean(post.diff^2,na.rm=T))
    
    ## Adding uncertainty estimates 
    
    if (se==1){
      
      XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
      
      for (k in 1:iterations){
        
        repeat{
          ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
          if (length(ctrls.k)>1){break}
        }
        
        dat.i.k=dataprep(D.i,
                         predictors=X1, 
                         predictors.op="mean",
                         dependent=Y1,
                         unit.variable="cid",
                         time.variable="year",
                         special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                                 list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                                 list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                         treatment.identifier=0,
                         controls.identifier=ctrls.k,
                         time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                         time.optimize.ssr=(max(1980,yr-6):yr),
                         unit.names.variable="Name",
                         time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
          dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
        }
        
        synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
        print(synth.tables)
        
        XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
      }
      
      if (yr<1986){
        synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
        synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
        synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
      } 
      
      if (yr>2008){
        synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
        synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
        synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
      } 
      
      if (yr>=1986 & yr<=2008){
        synthSE.TR[,1]<-dat.i.k$Y1plot
        synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
        synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
      }
      
      lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
      lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    }
    
    ## Figure A8a ##
    dev.off()
  
  
# Changes to the donor pool  
  
  ## Donor pool only IMF programs in the same year ## 
  SCM.R.1=function(treated, control, sunit, eunit, Y1, Y0, X1){
    
    for (i in sunit:eunit){
      cty=d$cid[treated[i]]
      yr=d$year[treated[i]]
      
      ctrls=intersect(control, which(d$spell==1 & d$year==yr))
      ctrls.cid=d$cid[ctrls]
      
      if (length(ctrls.cid)>1 & yr>1980){
        
        dat.i=dataprep(D,
                       predictors=X1,
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean")),
                       treatment.identifier=cty,
                       controls.identifier=ctrls.cid,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="iso3",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
          dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
        }
        
        synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        
        synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
        print(synth.tables)
        
        Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
        
        path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
        legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
        abline(v=yr, lty=2, col="gray")
        abline(h=0, lty=2, col="gray")
      }
    }
  }
  
  pdf("IDN-1997-R1.pdf")
  SCM.R.1(Treated.EXT, Control.EXT, idn, idn, "KOFEcGIdj", "trade_WDI", X1)
  ## Figure A8a ##
  dev.off()
  
  
  ## Donor pool including all observations but with IMF propensity score (imfPS) ##
  d$mingap.EXT=1
  MinGap=5
  for (i in 1:nrow(d)){
    idx.prior=which(d$iso3==d$iso3[i] & d$year<d$year[i] & d$year>=d$year[i]-MinGap)
    idx.after=which(d$iso3==d$iso3[i] & d$year>=d$year[i] & d$year<=d$year[i]+5)
    d$mingap.EXT[i]=as.numeric(all(d$IMFnn[idx.prior]==0))
  }
  
  SCM.R.2=function(treated, control, sunit, eunit, Y1, Y0, X1){
    
    for (i in sunit:eunit){
      cty=d$cid[treated[i]]
      yr=d$year[treated[i]]
      
      ctrls.cid=d$cid[d$mingap.EXT==1 & d$year==yr]
      ctrls.cid=ctrls.cid[is.na(ctrls.cid)==F]
      ctrls.cid=setdiff(ctrls.cid, cty)
      
      if (length(ctrls.cid)>1 & yr>1980){
        
        dat.i=dataprep(D,
                       predictors=X1,
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("imfPS",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean")),  
                       treatment.identifier=cty,
                       controls.identifier=ctrls.cid,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="iso3",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
          dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
        }
        
        synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        
        synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
        print(synth.tables)
        synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
        
        Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
        
        path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
        legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
        abline(v=yr, lty=2, col="gray")
        abline(h=0, lty=2, col="gray")
      }    
      
      ## add SEs ##
      XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
      
      for (k in 1:iterations){
        
        repeat{
          ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
          if (length(ctrls.k)>1){break}
        }
        
        dat.i.k=dataprep(D,
                       predictors=X1,
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("imfPS",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean")),  
                       treatment.identifier=cty,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="iso3",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
          dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
        }
        
        synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
        print(synth.tables)
        
        XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
      }
      
      if (yr<1986){
        #synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
        #synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
        synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
      } 
      
      if (yr>2008){
        #synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
        #synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
        synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
      } 
      
      if (yr>=1986 & yr<=2008){
        #synthSE.TR[,1]<-dat.i.k$Y1plot
        #synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
        synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
      }
      
      lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
      lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    }
  }
  
  X1P=c("imfPS",X1)
  pdf("IDN-1997-R2.pdf")
  SCM.R.2(Treated.EXT, Control.EXT, idn, idn, "KOFEcGIdj", "trade_WDI", X1P)
  ## Figure A8b ##
  dev.off()
  
  
  ## Donor pool additionally controlling for propensity score of sector conditionality (PS.EXT) ##
  i=idn
  treated=Treated.EXT 
  control=Control.EXT 
  Y1="KOFEcGIdj"
  Y0="trade_WDI"
  X1=c("PS.EXT", "lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  cty=d$cid[treated[i]]
  yr=1997
  D.i=CreateControlData(cty,yr,"Any")
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  
  pdf("IDN-1997-PSEXT.pdf")  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[1,1]<-d$iso3[treated[i]]
  synthSE.diag[1,2]<-yr
  synthSE.diag[1,3]<-pre.years 
  synthSE.diag[1,4]<-post.years 
  synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(
                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,1]<-dat.i.k$Y1plot
      synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
  }
  
  ## Figure A8c ##
  dev.off()
  
  
  ## Donor pool excluding countries from the same region ## 
  i=idn
  treated=Treated.EXT 
  control=Control.EXT 
  Y1="KOFEcGIdj"
  Y0="trade_WDI"
  X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  cty=d$cid[treated[i]]
  yr=1997
  D.i=CreateControlData(cty,yr,"Any")
  cid.reg=unique(D.i$cid[which(D.i$regid%in%D.i$regid[1])])
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  ctrls.cid.nreg=setdiff(ctrls.cid, cid.reg)
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid.nreg,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  
  pdf("IDN-1997-nreg.pdf")  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  legend("topright", paste("P=",length(ctrls.cid.nreg),sep=""), cex=0.55)
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[1,1]<-d$iso3[treated[i]]
  synthSE.diag[1,2]<-yr
  synthSE.diag[1,3]<-pre.years 
  synthSE.diag[1,4]<-post.years 
  synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid.nreg, length(ctrls.cid.nreg), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(
                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,1]<-dat.i.k$Y1plot
      synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
  }
  
  ## Figure A8d ## 
  dev.off()
  

  ## Donor pool of units with no program in the seven years prior to treatment ##
  MinGap=7
  Control.EXT7=Spells[1]
  Control.FIN7=numeric(0)
  
  for (s in 2:length(Spells)){
    diff=ifelse(d$cid[Spells[s]]==d$cid[Spells[s-1]], d$year[Spells[s]]-d$year[Spells[s-1]], 99)
    if (diff>=MinGap){
      Treated=c(Treated, Spells[s])
    }
    
    p.obs=ifelse(s<length(Spells), 
                 (Spells[s]:(Spells[s+1]-1))[d$IMFnn[Spells[s]:(Spells[s+1]-1)]==1],
                 Spells[s])
    
    if (any(d$SCsEXT[p.obs], na.rm=T)>0){
      #Treated.EXT7=c(Treated.EXT7, Spells[s])
    } else {
      Control.EXT7=c(Control.EXT7, Spells[s])
    }
    
    if (any(d$SCsFIN[p.obs], na.rm=T)>0){
      #Treated.FIN7=c(Treated.FIN7, Spells[s])
    } else {
      Control.FIN7=c(Control.FIN7, Spells[s])
    }
  }
  
  Control.EXT7=intersect(Treated, Control.EXT7)
  
  idn=16
  i=idn
  treated=Treated.EXT 
  control=Control.EXT7 
  Control.EXT=Control.EXT7
  Y1="KOFEcGIdj"
  Y0="trade_WDI"
  X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  cty=d$cid[treated[i]]
  yr=d$year[treated[i]]
  D.i=CreateControlData(cty,yr,"Any")
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  
  pdf("IDN-1997-R5.pdf")  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[1,1]<-d$iso3[treated[i]]
  synthSE.diag[1,2]<-yr
  synthSE.diag[1,3]<-pre.years 
  synthSE.diag[1,4]<-post.years 
  synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,1]<-dat.i.k$Y1plot
      synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
  }
  
  ## Figure A8e ##
  dev.off()
  
  
  ## Placebo treatment ten years later ## 
  cty=d$cid[treated[i]]
  yr=2007
  D.i=CreateControlData(cty,yr,"Any")
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  
  pdf("IDN-placebo2007.pdf")  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[1,1]<-d$iso3[treated[i]]
  synthSE.diag[1,2]<-yr
  synthSE.diag[1,3]<-pre.years 
  synthSE.diag[1,4]<-post.years 
  synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,1]<-dat.i.k$Y1plot
      synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
  }
  
  ## Figure A8f ##
  dev.off()
  
  
# Additional controls in the SCM matching 
  
  ## Political instability ##
  idn=16
  i=idn
  treated=Treated.EXT 
  control=Control.EXT 
  Y1="KOFEcGIdj"
  Y0="trade_WDI"
  X1=c("anylc_5y","lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  cty=d$cid[treated[i]]
  yr=d$year[treated[i]]
  D.i=CreateControlData(cty,yr,"Any")
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean"),
                                         list("lten",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  
  pdf("IDN-1997-polstab.pdf")  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[1,1]<-d$iso3[treated[i]]
  synthSE.diag[1,2]<-yr
  synthSE.diag[1,3]<-pre.years 
  synthSE.diag[1,4]<-post.years 
  synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,1]<-dat.i.k$Y1plot
      synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
  }
  
  ## Figure A9a ##
  dev.off()
  
  ## Business cycle effects ##
  avgdp=as.data.frame(D%>% group_by(iso3) %>% summarise(avgdp=mean(gdp_growth_WDI)))
  D$bcycle=NA
  for (k in 1:nrow(avgdp)){
    idx=which(D$iso3%in%avgdp$iso3[k])
    D$bcycle[idx]=D$gdp_growth_WDI[idx]-avgdp$avgdp[k]
  }
  
  X1=c("bcycle","lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  cty=d$cid[treated[i]]
  yr=d$year[treated[i]]
  D.i=CreateControlData(cty,yr,"Any")
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean"),
                                         list("lten",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  
  pdf("IDN-1997-bcycle.pdf")  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[1,1]<-d$iso3[treated[i]]
  synthSE.diag[1,2]<-yr
  synthSE.diag[1,3]<-pre.years 
  synthSE.diag[1,4]<-post.years 
  synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,1]<-dat.i.k$Y1plot
      synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
  }
  
  ## Figure A9b ## 
  dev.off() 
  
  
  ## IMF Technical Assistance eligibility ##
  X1=c("imf_ro","lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  cty=d$cid[treated[i]]
  yr=d$year[treated[i]]
  D.i=CreateControlData(cty,yr,"Any")
  
  ctrls.cid=setdiff(unique(D.i$cid), 0)
  
  dat.i=dataprep(D.i,
                 predictors=X1, 
                 predictors.op="mean",
                 dependent=Y1,
                 unit.variable="cid",
                 time.variable="year",
                 special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                         list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                         list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                 treatment.identifier=0,
                 controls.identifier=ctrls.cid,
                 time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                 time.optimize.ssr=(max(1980,yr-6):yr),
                 unit.names.variable="Name",
                 time.plot=(max(1980,yr-6):min(2014,yr+6))
  )
  
  if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
    idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
    dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
  }
  
  synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
  if (inherits(synth.out, "error")) next
  
  synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
  print(synth.tables)
  synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
  synth.tables$tab.v
  
  Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
  
  pdf("IDN-1997-TA.pdf")  
  path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
  legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
  abline(v=yr, lty=2, col="gray")
  abline(h=0, lty=2, col="gray")
  
  pre.years=6+(min(yr,1986)-1986)
  post.years=6-(max(yr,2008)-2008)
  post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
  pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
  
  synthSE.diag[1,1]<-d$iso3[treated[i]]
  synthSE.diag[1,2]<-yr
  synthSE.diag[1,3]<-pre.years 
  synthSE.diag[1,4]<-post.years 
  synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
  synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
  
  if (se==1){
    
    XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
    
    for (k in 1:iterations){
      
      repeat{
        ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
        if (length(ctrls.k)>1){break}
      }
      
      dat.i.k=dataprep(D.i,
                       predictors=X1, 
                       predictors.op="mean",
                       dependent=Y1,
                       unit.variable="cid",
                       time.variable="year",
                       special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                               list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                               list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                       treatment.identifier=0,
                       controls.identifier=ctrls.k,
                       time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                       time.optimize.ssr=(max(1980,yr-6):yr),
                       unit.names.variable="Name",
                       time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
        dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
      print(synth.tables)
      
      XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
    }
    
    if (yr<1986){
      synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
      synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
      synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
    } 
    
    if (yr>2008){
      synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
      synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
      synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
    } 
    
    if (yr>=1986 & yr<=2008){
      synthSE.TR[,1]<-dat.i.k$Y1plot
      synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
      synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
    }
    
    lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
  }
  
  ## Figure A9c ## 
  dev.off() 
  
  X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
  
  
  
# Robustness checks 
###################################################################################
  
# SCM routine with 'political instability' as additional  control 
  
  ## EXT sector ##
  
  pdf("EXT-EcGI-SE-polstab.pdf")
  
  synthSE.TR=matrix(NA, ncol=length(Treated.EXT), nrow=13)
  synthSE.CR=synthSE.TR
  synthSE.se=synthSE.TR
  synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
  synthSE.diff=expand.grid(iso3=d$iso3[Treated.EXT], t=0:6, diff=NA)
  synthSE.weights.var=matrix(NA, ncol=length(Treated.EXT), nrow=18+2)
  
  for (i in 1:length(Treated.EXT)){
    
    treated=Treated.EXT 
    control=Control.EXT 
    Y1="KOFEcGIdj"
    Y0="trade_WDI"
    X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
    
    cty=d$cid[treated[i]]
    yr=d$year[treated[i]]
    D.i=CreateControlData(cty,yr,"Any")
    
    ## Matching
    
    ctrls.cid=setdiff(unique(D.i$cid), 0)
    
    ## Main routine  
    
    dat.i=dataprep(D.i,
                   predictors=X1, 
                   predictors.op="mean",
                   dependent=Y1,
                   unit.variable="cid",
                   time.variable="year",
                   special.predictors=list(list("anylc_10y5y",yr,"mean"),list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                           list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                           list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                   treatment.identifier=0,
                   controls.identifier=ctrls.cid,
                   time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                   time.optimize.ssr=(max(1980,yr-6):yr),
                   unit.names.variable="Name",
                   time.plot=(max(1980,yr-6):min(2014,yr+6))
    )
    
    if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
      idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
      dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
    }
    
    synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
    if (inherits(synth.out, "error")) next
    
    synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
    print(synth.tables)
    synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
    
    Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
    
    path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
    legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
    abline(v=yr, lty=2, col="gray")
    abline(h=0, lty=2, col="gray")
    
    pre.years=6+(min(yr,1986)-1986)
    post.years=6-(max(yr,2008)-2008)
    post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
    pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
    
    synthSE.diag[1,1]<-d$iso3[treated[i]]
    synthSE.diag[1,2]<-yr
    synthSE.diag[1,3]<-pre.years 
    synthSE.diag[1,4]<-post.years 
    synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
    synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
    synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
    
    if (se==1){
      
      XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
      
      for (k in 1:iterations){
        
        repeat{
          ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
          if (length(ctrls.k)>1){break}
        }
        
        dat.i.k=dataprep(D.i,
                         predictors=X1, 
                         predictors.op="mean",
                         dependent=Y1,
                         unit.variable="cid",
                         time.variable="year",
                         special.predictors=list(list("anylc_10y5y",yr,"mean"),list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                                 list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                                 list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                         treatment.identifier=0,
                         controls.identifier=ctrls.k,
                         time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                         time.optimize.ssr=(max(1980,yr-6):yr),
                         unit.names.variable="Name",
                         time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
          dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
        }
        
        synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
        print(synth.tables)
        
        XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
      }
      
      if (yr<1986){
        synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
        synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
        synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
      } 
      
      if (yr>2008){
        synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
        synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
        synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
      } 
      
      if (yr>=1986 & yr<=2008){
        synthSE.TR[,1]<-dat.i.k$Y1plot
        synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
        synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
      }
      
      lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
      lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    }
  }
    
  write.csv(synthSE.diag, "EXT-EcGI-SE_1-polstab.csv")
  write.csv(synthSE.TR, "EXT-EcGI-SE_2-polstab.csv")
  write.csv(synthSE.CR, "EXT-EcGI-SE_3-polstab.csv")
  write.csv(synthSE.se, "EXT-EcGI-SE_4-polstab.csv")  
  write.csv(synthSE.diff, "EXT-EcGI-SE_5-polstab.csv")
  write.csv(synthSE.weights.var, "EXT-EcGI-w-EXT-polstab.csv")
    
  save.image(file="R-EXT-EcGI-polstab.Rdata")
    
    
    ## FIN sector ## 
    
    synthSE.TR=matrix(NA, ncol=length(Treated.FIN), nrow=13)
    synthSE.CR=synthSE.TR
    synthSE.se=synthSE.TR
    synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
    synthSE.diff=expand.grid(iso3=d$iso3[Treated.FIN], t=0:6, diff=NA)
    synthSE.weights.var=matrix(NA, ncol=length(Treated.FIN), nrow=18+2)
    
    pdf("FIN-FiGI-SE-polstab.pdf")
    
    for (i in 1:length(Treated.FIN)){
      
      treated=Treated.FIN
      control=Control.FIN 
      Y1="KOFFiGIdj"
      Y0="fdi_gdp_WDI"
      X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
      
      cty=d$cid[treated[i]]
      yr=d$year[treated[i]]
      D.i=CreateControlData(cty,yr,"Any")
      
      ## Matching
      
      ctrls.cid=setdiff(unique(D.i$cid), 0)
      
      ## Main routine  
      
      dat.i=dataprep(D.i,
                     predictors=X1, 
                     predictors.op="mean",
                     dependent=Y1,
                     unit.variable="cid",
                     time.variable="year",
                     special.predictors=list(list("anylc_10y5y",yr,"mean"),list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                             list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                             list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                     treatment.identifier=0,
                     controls.identifier=ctrls.cid,
                     time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                     time.optimize.ssr=(max(1980,yr-6):yr),
                     unit.names.variable="Name",
                     time.plot=(max(1980,yr-6):min(2014,yr+6))
      )
      
      if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
        idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
        dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
      }
      
      synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
      
      synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
      print(synth.tables)
      
      Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
      path.plot(synth.res=synth.out, dataprep.res=dat.i, Main=paste(d$cname[treated[i]],yr), Ylab="De-jure financial liberalization", Xlab=paste("|D|=",length(ctrls.cid),sep=""), 
                Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
      abline(v=yr, lty=2, col="gray")
      abline(h=0, lty=2, col="gray")
      
      pre.years=6+(min(yr,1986)-1986)
      post.years=6-(max(yr,2008)-2008)
      post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
      pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
      
      synthSE.diag[i,1]<-d$iso3[treated[i]]
      synthSE.diag[i,2]<-yr
      synthSE.diag[i,3]<-pre.years 
      synthSE.diag[i,4]<-post.years 
      synthSE.diag[i,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
      synthSE.diag[i,6]<-sqrt(mean(pre.diff^2,na.rm=T))
      synthSE.diag[i,7]<-sqrt(mean(post.diff^2,na.rm=T))
      
      synthSE.diff$diff[(synthSE.diff$iso3 %in% d$iso3[treated[i]])[1:(post.years+1)]]<-post.diff
      synthSE.weights.var[1:19,i]=as.numeric(synth.tables$tab.v)
      synthSE.weights.var[20,i]=mahalanobis(synth.tables$tab.pred[1:12,1],synth.tables$tab.pred[1:12,2],cov=cov(D[,X1]))
      
      ## Adding uncertainty estimates 
      
      if (se==1){
        
        XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
        
        for (k in 1:iterations){
          
          repeat{
            ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
            if (length(ctrls.k)>1){break}
          }
          
          dat.i.k=dataprep(D.i,
                           predictors=X1, 
                           predictors.op="mean",
                           dependent=Y1,
                           unit.variable="cid",
                           time.variable="year",
                           special.predictors=list(list("anylc_10y5y",yr,"mean"),list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                                   list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                                   list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                           treatment.identifier=0,
                           controls.identifier=ctrls.k,
                           time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                           time.optimize.ssr=(max(1980,yr-6):yr),
                           unit.names.variable="Name",
                           time.plot=(max(1980,yr-6):min(2014,yr+6))
          )
          
          if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
            idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
            dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
          }
          
          synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
          if (inherits(synth.out, "error")) next
          synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
          print(synth.tables)
          
          XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
        }
        
        if (yr<1986){
          synthSE.TR[,i]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
          synthSE.CR[,i]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
          synthSE.se[,i]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
        } 
        
        if (yr>2008){
          synthSE.TR[,i]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
          synthSE.CR[,i]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
          synthSE.se[,i]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
        } 
        
        if (yr>=1986 & yr<=2008){
          synthSE.TR[,i]<-dat.i.k$Y1plot
          synthSE.CR[,i]<-apply(XR,1,mean,na.rm=T)
          synthSE.se[,i]<-sqrt(apply(XR, 1, var, na.rm=T))
        }
        
        lines(rownames(dat.i$Y1plot), Y0.w.synth+1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
        lines(rownames(dat.i$Y1plot), Y0.w.synth-1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
      }
    }
    
    dev.off()
    write.csv(synthSE.diag, "FIN-FiGI-SE_1-polstab.csv")
    write.csv(synthSE.TR, "FIN-FiGI-SE_2-polstab.csv")
    write.csv(synthSE.CR, "FIN-FiGI-SE_3-polstab.csv")
    write.csv(synthSE.se, "FIN-FiGI-SE_4-polstab.csv")  
    write.csv(synthSE.diff, "FIN-FiGI-SE_5-polstab.csv")
    write.csv(synthSE.weights.var, "FIN-FiGI-w-polstab.csv")
    
    save.image(file="R-FIN-FiGI-polstab.Rdata")
    
    
# Three-year gap rule   
    
  MinGap=3
  Control.EXT3=Spells[1]
  Control.FIN3=numeric(0)
  Treated.EXT3=numeric(0)
  Treated.FIN3=Spells[1]
  
  for (s in 2:length(Spells)){
    diff=ifelse(d$cid[Spells[s]]==d$cid[Spells[s-1]], d$year[Spells[s]]-d$year[Spells[s-1]], 99)
    if (diff>=MinGap){
      Treated=c(Treated, Spells[s])
    }
    
    p.obs=ifelse(s<length(Spells), 
                 (Spells[s]:(Spells[s+1]-1))[d$IMFnn[Spells[s]:(Spells[s+1]-1)]==1],
                 Spells[s])
    
    if (any(d$SCsEXT[p.obs], na.rm=T)>0){
      Treated.EXT3=c(Treated.EXT3, Spells[s])
      } else {
      Control.EXT3=c(Control.EXT3, Spells[s])
      }
    
    if (any(d$SCsFIN[p.obs], na.rm=T)>0){
      Treated.FIN3=c(Treated.FIN3, Spells[s])
      } else {
      Control.FIN3=c(Control.FIN3, Spells[s])
      }
    }
  
  Treated.EXT3=intersect(Treated, Treated.EXT3)
  Control.EXT3=intersect(Treated, Control.EXT3)
  Treated.FIN3=intersect(Treated, Treated.FIN3)
  Control.FIN3=intersect(Treated, Control.FIN3)  
    
  ## EXT sector ## 
  synthSE.TR=matrix(NA, ncol=length(Treated.EXT3), nrow=13)
  synthSE.CR=synthSE.TR
  synthSE.se=synthSE.TR
  synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
  synthSE.diff=expand.grid(iso3=d$iso3[Treated.EXT3], t=0:6, diff=NA)
  synthSE.weights.var=matrix(NA, ncol=length(Treated.EXT3), nrow=18+1)
  
  pdf("EXT3-EcGI-SE.pdf")
    
  for (i in 1:length(Treated.EXT3)){
      
    treated=Treated.EXT3 
    control=Control.EXT3 
    Y1="KOFEcGIdj"
    Y0="trade_WDI"
    X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y")
    
    cty=d$cid[treated[i]]
    yr=d$year[treated[i]]
    D.i=CreateControlData(cty,yr,"Any")
      
    ## Matching
      
    ctrls.cid=setdiff(unique(D.i$cid), 0)
    
    ## Main routine  
    
    dat.i=dataprep(D.i,
                   predictors=X1, 
                   predictors.op="mean",
                   dependent=Y1,
                   unit.variable="cid",
                   time.variable="year",
                   special.predictors=list(list("anylc_10y5y",yr,"mean"),list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                           list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                           list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                   treatment.identifier=0,
                   controls.identifier=ctrls.cid,
                   time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                   time.optimize.ssr=(max(1980,yr-6):yr),
                   unit.names.variable="Name",
                   time.plot=(max(1980,yr-6):min(2014,yr+6))
    )
    
    if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
      idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
      dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
    }
    
    synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
    if (inherits(synth.out, "error")) next
    
    synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
    print(synth.tables)
    
    Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
    path.plot(synth.res=synth.out, dataprep.res=dat.i, Main=paste(d$cname[treated[i]],yr), Ylab="De-jure economic liberalization", Xlab=paste("|D|=",length(ctrls.cid),sep=""), 
              Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
    abline(v=yr, lty=2, col="gray")
    abline(h=0, lty=2, col="gray")
    
    pre.years=6+(min(yr,1986)-1986)
    post.years=6-(max(yr,2008)-2008)
    post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
    pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
    
    synthSE.diag[i,1]<-d$iso3[treated[i]]
    synthSE.diag[i,2]<-yr
    synthSE.diag[i,3]<-pre.years 
    synthSE.diag[i,4]<-post.years 
    synthSE.diag[i,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
    synthSE.diag[i,6]<-sqrt(mean(pre.diff^2,na.rm=T))
    synthSE.diag[i,7]<-sqrt(mean(post.diff^2,na.rm=T))
    
    synthSE.diff$diff[(synthSE.diff$iso3 %in% d$iso3[treated[i]])[1:(post.years+1)]]<-post.diff
    synthSE.weights.var[1:18,i]=as.numeric(synth.tables$tab.v)
    synthSE.weights.var[19,i]=mahalanobis(synth.tables$tab.pred[1:11,1],synth.tables$tab.pred[1:11,2],cov=cov(D[,X1]))
    
    ## Adding uncertainty estimates 
    
    if (se==1){
      
      XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
      
      for (k in 1:iterations){
        
        repeat{
          ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
          if (length(ctrls.k)>1){break}
        }
        
        dat.i.k=dataprep(D.i,
                         predictors=X1, 
                         predictors.op="mean",
                         dependent=Y1,
                         unit.variable="cid",
                         time.variable="year",
                         special.predictors=list(list("anylc_10y5y",yr,"mean"),list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                                 list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                                 list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                         treatment.identifier=0,
                         controls.identifier=ctrls.k,
                         time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                         time.optimize.ssr=(max(1980,yr-6):yr),
                         unit.names.variable="Name",
                         time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
          dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
        }
        
      synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
        print(synth.tables)
        
        XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
      }
      
      if (yr<1986){
        synthSE.TR[,i]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
        synthSE.CR[,i]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
        synthSE.se[,i]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
      } 
      
      if (yr>2008){
        synthSE.TR[,i]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
        synthSE.CR[,i]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
        synthSE.se[,i]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
      } 
      
      if (yr>=1986 & yr<=2008){
        synthSE.TR[,i]<-dat.i.k$Y1plot
        synthSE.CR[,i]<-apply(XR,1,mean,na.rm=T)
        synthSE.se[,i]<-sqrt(apply(XR, 1, var, na.rm=T))
      }
      
      lines(rownames(dat.i$Y1plot), Y0.w.synth+1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
      lines(rownames(dat.i$Y1plot), Y0.w.synth-1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
    }
  }
    
  dev.off()
  write.csv(synthSE.diag, "EXT3-EcGI-SE_1.csv")
  write.csv(synthSE.TR, "EXT3-EcGI-SE_2.csv")
  write.csv(synthSE.CR, "EXT3-EcGI-SE_3.csv")
  write.csv(synthSE.se, "EXT3-EcGI-SE_4.csv")  
  write.csv(synthSE.diff, "EXT3-EcGI-SE_5.csv")
  write.csv(synthSE.weights.var, "EXT3-EcGI-w-EXT.csv")
  
  save.image(file="R-EXT3-EcGI.Rdata")
    
  ## FIN ##  
  synthSE.TR=matrix(NA, ncol=length(Treated.FIN3), nrow=13)
  synthSE.CR=synthSE.TR
  synthSE.se=synthSE.TR
  synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
  synthSE.diff=expand.grid(iso3=d$iso3[Treated.FIN3], t=0:6, diff=NA)
  synthSE.weights.var=matrix(NA, ncol=length(Treated.FIN3), nrow=18+1)
  
  pdf("FIN3-EcGI-SE.pdf")
  
  for (i in 1:length(Treated.FIN3)){
    
    treated=Treated.FIN3
    control=Control.FIN3
    #Y1="KOFFiGIdj"
    Y0="fdi_gdp_WDI"
    X1=c("lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y")
    
    cty=d$cid[treated[i]]
    yr=d$year[treated[i]]
    D.i=CreateControlData(cty,yr,"Any")
    
    ## Matching
      
    ctrls.cid=setdiff(unique(D.i$cid), 0)
    
    ## Main routine  
      
    dat.i=dataprep(D.i,
                   predictors=X1, 
                   predictors.op="mean",
                   dependent=Y1,
                   unit.variable="cid",
                   time.variable="year",
                   special.predictors=list(list("anylc_10y5y",yr,"mean"),list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                           list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                           list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                   treatment.identifier=0,
                   controls.identifier=ctrls.cid,
                   time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                   time.optimize.ssr=(max(1980,yr-6):yr),
                   unit.names.variable="Name",
                   time.plot=(max(1980,yr-6):min(2014,yr+6))
    )
    
    if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
      idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
      dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
    }
    
    synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
      if (inherits(synth.out, "error")) next
  
    synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
    print(synth.tables)
    
    Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
    path.plot(synth.res=synth.out, dataprep.res=dat.i, Main=paste(d$cname[treated[i]],yr), Ylab="De-jure economic liberalization", Xlab=paste("|D|=",length(ctrls.cid),sep=""), 
              Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
    abline(v=yr, lty=2, col="gray")
    abline(h=0, lty=2, col="gray")
    
    pre.years=6+(min(yr,1986)-1986)
    post.years=6-(max(yr,2008)-2008)
    post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
    pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
    
    synthSE.diag[i,1]<-d$iso3[treated[i]]
    synthSE.diag[i,2]<-yr
    synthSE.diag[i,3]<-pre.years 
    synthSE.diag[i,4]<-post.years 
    synthSE.diag[i,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
    synthSE.diag[i,6]<-sqrt(mean(pre.diff^2,na.rm=T))
    synthSE.diag[i,7]<-sqrt(mean(post.diff^2,na.rm=T))
    
    synthSE.diff$diff[(synthSE.diff$iso3 %in% d$iso3[treated[i]])[1:(post.years+1)]]<-post.diff
    synthSE.weights.var[1:18,i]=as.numeric(synth.tables$tab.v)
    synthSE.weights.var[19,i]=mahalanobis(synth.tables$tab.pred[1:11,1],synth.tables$tab.pred[1:11,2],cov=cov(D[,X1]))
    
    ## Adding uncertainty estimates 
    
    if (se==1){
      
      XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
      
      for (k in 1:iterations){
        
        repeat{
          ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
          if (length(ctrls.k)>1){break}
        }
        
        dat.i.k=dataprep(D.i,
                         predictors=X1, 
                         predictors.op="mean",
                         dependent=Y1,
                         unit.variable="cid",
                         time.variable="year",
                         special.predictors=list(list("anylc_10y5y",yr,"mean"),list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                                 list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                                 list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                         treatment.identifier=0,
                         controls.identifier=ctrls.k,
                         time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                         time.optimize.ssr=(max(1980,yr-6):yr),
                         unit.names.variable="Name",
                         time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
          dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
        }
        
        synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
        print(synth.tables)
        
        XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
      }
      
      if (yr<1986){
        synthSE.TR[,i]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
        synthSE.CR[,i]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
        synthSE.se[,i]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
      } 
      
      if (yr>2008){
        synthSE.TR[,i]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
        synthSE.CR[,i]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
        synthSE.se[,i]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
      } 
      
      if (yr>=1986 & yr<=2008){
        synthSE.TR[,i]<-dat.i.k$Y1plot
        synthSE.CR[,i]<-apply(XR,1,mean,na.rm=T)
        synthSE.se[,i]<-sqrt(apply(XR, 1, var, na.rm=T))
      }
      
      lines(rownames(dat.i$Y1plot), Y0.w.synth+1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
      lines(rownames(dat.i$Y1plot), Y0.w.synth-1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
    }
  }
  
  dev.off()
  write.csv(synthSE.diag, "FIN3-FiGI-SE_1.csv")
  write.csv(synthSE.TR, "FIN3-FiGI-SE_2.csv")
  write.csv(synthSE.CR, "FIN3-FiGI-SE_3.csv")
  write.csv(synthSE.se, "FIN3-FiGI-SE_4.csv")  
  write.csv(synthSE.diff, "FIN3-FiGI-SE_5.csv")
  write.csv(synthSE.weights.var, "FIN3-FiGI-w.csv")
  
  save.image(file="R-FIN3-FiGI.Rdata")
    
    
# Controlling for IMF Technical Assistance   
    
  ## EXT sector ##
  
  pdf("EXT-EcGI-SE-TA.pdf")
  
  synthSE.TR=matrix(NA, ncol=length(Treated.EXT), nrow=13)
  synthSE.CR=synthSE.TR
  synthSE.se=synthSE.TR
  synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
  synthSE.diff=expand.grid(iso3=d$iso3[Treated.EXT], t=0:6, diff=NA)
  synthSE.weights.var=matrix(NA, ncol=length(Treated.EXT), nrow=18+2)
  
  for (i in 1:length(Treated.EXT)){
    
    treated=Treated.EXT 
    control=Control.EXT 
    Y1="KOFEcGIdj"
    Y0="trade_WDI"
    X1=c("imf_ro","lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
    
    cty=d$cid[treated[i]]
    yr=d$year[treated[i]]
    D.i=CreateControlData(cty,yr,"Any")
    
    ## Matching
    
    ctrls.cid=setdiff(unique(D.i$cid), 0)
    
    ## Main routine  
    
    dat.i=dataprep(D.i,
                   predictors=X1, 
                   predictors.op="mean",
                   dependent=Y1,
                   unit.variable="cid",
                   time.variable="year",
                   special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                           list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                           list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                   treatment.identifier=0,
                   controls.identifier=ctrls.cid,
                   time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                   time.optimize.ssr=(max(1980,yr-6):yr),
                   unit.names.variable="Name",
                   time.plot=(max(1980,yr-6):min(2014,yr+6))
    )
    
    if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
      idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
      dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
    }
    
    synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
    if (inherits(synth.out, "error")) next
    
    synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
    print(synth.tables)
    synth.tables$tab.w[synth.tables$tab.w$w.weights>0,]
    
    Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
    
    path.plot(synth.res=synth.out, dataprep.res=dat.i, Ylab="De-jure economic liberalization", Xlab="", Main=NA, Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
    legend("topright", paste("P=",length(ctrls.cid),sep=""), cex=0.55)
    abline(v=yr, lty=2, col="gray")
    abline(h=0, lty=2, col="gray")
    
    pre.years=6+(min(yr,1986)-1986)
    post.years=6-(max(yr,2008)-2008)
    post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
    pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
    
    synthSE.diag[1,1]<-d$iso3[treated[i]]
    synthSE.diag[1,2]<-yr
    synthSE.diag[1,3]<-pre.years 
    synthSE.diag[1,4]<-post.years 
    synthSE.diag[1,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
    synthSE.diag[1,6]<-sqrt(mean(pre.diff^2,na.rm=T))
    synthSE.diag[1,7]<-sqrt(mean(post.diff^2,na.rm=T))
    
    if (se==1){
      
      XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
      
      for (k in 1:iterations){
        
        repeat{
          ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
          if (length(ctrls.k)>1){break}
        }
        
        dat.i.k=dataprep(D.i,
                         predictors=X1, 
                         predictors.op="mean",
                         dependent=Y1,
                         unit.variable="cid",
                         time.variable="year",
                         special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                                 list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                                 list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                         treatment.identifier=0,
                         controls.identifier=ctrls.k,
                         time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                         time.optimize.ssr=(max(1980,yr-6):yr),
                         unit.names.variable="Name",
                         time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
          dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
        }
        
        synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
        print(synth.tables)
        
        XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
      }
      
      if (yr<1986){
        synthSE.TR[,1]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
        synthSE.CR[,1]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
        synthSE.se[,1]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
      } 
      
      if (yr>2008){
        synthSE.TR[,1]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
        synthSE.CR[,1]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
        synthSE.se[,1]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
      } 
      
      if (yr>=1986 & yr<=2008){
        synthSE.TR[,1]<-dat.i.k$Y1plot
        synthSE.CR[,1]<-apply(XR,1,mean,na.rm=T)
        synthSE.se[,1]<-sqrt(apply(XR, 1, var, na.rm=T))
      }
      
      lines(rownames(dat.i$Y1plot), Y0.w.synth+1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
      lines(rownames(dat.i$Y1plot), Y0.w.synth-1.97*synthSE.se[is.na(synthSE.se[,1])==F, 1], lwd=1, lty=2)
    }
  }
  
  write.csv(synthSE.diag, "EXT-EcGI-SE_1-TA.csv")
  write.csv(synthSE.TR, "EXT-EcGI-SE_2-TA.csv")
  write.csv(synthSE.CR, "EXT-EcGI-SE_3-TA.csv")
  write.csv(synthSE.se, "EXT-EcGI-SE_4-TA.csv")  
  write.csv(synthSE.diff, "EXT-EcGI-SE_5-TA.csv")
  write.csv(synthSE.weights.var, "EXT-EcGI-w-EXT-TA.csv")
  
  save.image(file="R-EXT-EcGI-TA.Rdata")
  
  
  ## FIN sector ## 
  
  synthSE.TR=matrix(NA, ncol=length(Treated.FIN), nrow=13)
  synthSE.CR=synthSE.TR
  synthSE.se=synthSE.TR
  synthSE.diag=data.frame(iso3=NA, year=NA, pre.years=NA, post.years=NA, rmspe=NA, rmspe0=NA, rmspe1=NA)
  synthSE.diff=expand.grid(iso3=d$iso3[Treated.FIN], t=0:6, diff=NA)
  synthSE.weights.var=matrix(NA, ncol=length(Treated.FIN), nrow=18+2)
  
  pdf("FIN-FiGI-SE-TA.pdf")
  
  for (i in 2:length(Treated.FIN)){
    
    treated=Treated.FIN
    control=Control.FIN 
    Y1="KOFFiGIdj"
    Y0="fdi_gdp_WDI"
    X1=c("imf_ro","lngdppc","lnpop","capacity","globpol_KOF","gdp_growth_WDI","reserves_WDI","cab_gdp_WDI","debtser_gni_WDI","fuel_exp_WDI","polconiii","anym_5y","milexpgdp")
    
    cty=d$cid[treated[i]]
    yr=d$year[treated[i]]
    D.i=CreateControlData(cty,yr,"Any")
    
    ## Matching
    
    ctrls.cid=setdiff(unique(D.i$cid), 0)
    
    ## Main routine  
    
    dat.i=dataprep(D.i,
                   predictors=X1, 
                   predictors.op="mean",
                   dependent=Y1,
                   unit.variable="cid",
                   time.variable="year",
                   special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                           list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                           list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                   treatment.identifier=0,
                   controls.identifier=ctrls.cid,
                   time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                   time.optimize.ssr=(max(1980,yr-6):yr),
                   unit.names.variable="Name",
                   time.plot=(max(1980,yr-6):min(2014,yr+6))
    )
    
    if (any(apply(dat.i$X0, 1, var, na.rm=T)==0)==T){
      idx=which(apply(dat.i$X0, 1, var, na.rm=T)==0)
      dat.i$X0[idx,1]=dat.i$X0[idx,1]+runif(1)
    }
    
    synth.out=tryCatch(synth(dat.i, Margin.ipop=7e-3), error=function(e) e)
    if (inherits(synth.out, "error")) next
    
    synth.tables=synth.tab(dataprep.res=dat.i, synth.res=synth.out)
    print(synth.tables)
    
    Y0.w.synth=dat.i$Y0plot %*% synth.out$solution.w  
    path.plot(synth.res=synth.out, dataprep.res=dat.i, Main=paste(d$cname[treated[i]],yr), Ylab="De-jure financial liberalization", Xlab=paste("|D|=",length(ctrls.cid),sep=""), 
              Legend.position="topleft", Ylim=c(Lim(Y0.w.synth, dat.i$Y1plot, 0), Lim(Y0.w.synth, dat.i$Y1plot, 1)))
    abline(v=yr, lty=2, col="gray")
    abline(h=0, lty=2, col="gray")
    
    pre.years=6+(min(yr,1986)-1986)
    post.years=6-(max(yr,2008)-2008)
    post.diff=(dat.i$Y1plot-Y0.w.synth)[(length(dat.i$Y1plot)-post.years):(length(dat.i$Y1plot))]
    pre.diff=(dat.i$Y1plot-Y0.w.synth)[1:pre.years]
    
    synthSE.diag[i,1]<-d$iso3[treated[i]]
    synthSE.diag[i,2]<-yr
    synthSE.diag[i,3]<-pre.years 
    synthSE.diag[i,4]<-post.years 
    synthSE.diag[i,5]<-sqrt(mean(post.diff^2,na.rm=T)/mean(pre.diff^2,na.rm=T))
    synthSE.diag[i,6]<-sqrt(mean(pre.diff^2,na.rm=T))
    synthSE.diag[i,7]<-sqrt(mean(post.diff^2,na.rm=T))
    
    synthSE.diff$diff[(synthSE.diff$iso3 %in% d$iso3[treated[i]])[1:(post.years+1)]]<-post.diff
    synthSE.weights.var[1:19,i]=as.numeric(synth.tables$tab.v)
    synthSE.weights.var[20,i]=mahalanobis(synth.tables$tab.pred[1:12,1],synth.tables$tab.pred[1:12,2],cov=cov(D[,X1]))
    
    ## Adding uncertainty estimates 
    
    if (se==1){
      
      XR=matrix(NA, nrow=13-(max(yr,2008)-2008)+(min(yr,1986)-1986), ncol=iterations)
      
      for (k in 1:iterations){
        
        repeat{
          ctrls.k=sort(unique(sample(ctrls.cid, length(ctrls.cid), replace=T)))
          if (length(ctrls.k)>1){break}
        }
        
        dat.i.k=dataprep(D.i,
                         predictors=X1, 
                         predictors.op="mean",
                         dependent=Y1,
                         unit.variable="cid",
                         time.variable="year",
                         special.predictors=list(list("BA2TOT", yr,"mean"),list("BA2scope",yr,"mean"),
                                                 list(Y1,max(1980,yr-1),"mean"),list(Y1,max(1980,yr-5),"mean"),
                                                 list(Y0,max(1980,yr-1),"mean"),list("nUnder",max(1980,yr-1),"mean")),
                         treatment.identifier=0,
                         controls.identifier=ctrls.k,
                         time.predictors.prior=(max(1980,yr-6):max(1980,yr-1)),
                         time.optimize.ssr=(max(1980,yr-6):yr),
                         unit.names.variable="Name",
                         time.plot=(max(1980,yr-6):min(2014,yr+6))
        )
        
        if (any(apply(dat.i.k$X0, 1, var, na.rm=T)==0)==T){
          idx=which(apply(dat.i.k$X0, 1, var, na.rm=T)==0)
          dat.i.k$X0[idx,1]=dat.i.k$X0[idx,1]+runif(1)
        }
        
        synth.out=tryCatch(synth(dat.i.k, Margin.ipop=7e-3), error=function(e) e)
        if (inherits(synth.out, "error")) next
        synth.tables=synth.tab(dataprep.res=dat.i.k, synth.res=synth.out)
        print(synth.tables)
        
        XR[,k]=dat.i.k$Y0plot %*% synth.out$solution.w
      }
      
      if (yr<1986){
        synthSE.TR[,i]<-c(rep(NA,1986-yr), dat.i.k$Y1plot)
        synthSE.CR[,i]<-c(rep(NA,1986-yr), apply(XR,1,mean,na.rm=T))
        synthSE.se[,i]<-c(rep(NA,1986-yr), sqrt(apply(XR, 1, var, na.rm=T)))
      } 
      
      if (yr>2008){
        synthSE.TR[,i]<-c(dat.i.k$Y1plot, rep(NA,yr-2008))
        synthSE.CR[,i]<-c(apply(XR,1,mean,na.rm=T), rep(NA,yr-2008))
        synthSE.se[,i]<-c(sqrt(apply(XR, 1, var, na.rm=T)), rep(NA,yr-2008))
      } 
      
      if (yr>=1986 & yr<=2008){
        synthSE.TR[,i]<-dat.i.k$Y1plot
        synthSE.CR[,i]<-apply(XR,1,mean,na.rm=T)
        synthSE.se[,i]<-sqrt(apply(XR, 1, var, na.rm=T))
      }
      
      lines(rownames(dat.i$Y1plot), Y0.w.synth+1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
      lines(rownames(dat.i$Y1plot), Y0.w.synth-1.645*synthSE.se[is.na(synthSE.se[,i])==F, i], lwd=1, lty=2)
    }
  }
  
  dev.off()
  write.csv(synthSE.diag, "FIN-FiGI-SE_1-TA.csv")
  write.csv(synthSE.TR, "FIN-FiGI-SE_2-TA.csv")
  write.csv(synthSE.CR, "FIN-FiGI-SE_3-TA.csv")
  write.csv(synthSE.se, "FIN-FiGI-SE_4-TA.csv")  
  write.csv(synthSE.diff, "FIN-FiGI-SE_5-TA.csv")
  write.csv(synthSE.weights.var, "FIN-FiGI-w-TA.csv")
  
  
  #
#############
  #NOT USED#
###########
  
  setwd("C:\\berni\\forschung\\imfown\\data")
  
  d=read.dta("IMFOWN_00.dta")
  
  d=d[,c("cname","cid","iso3","year","region_WDI","KOFTrGIdj","KOFTrGIdf","KOFFiGIdj","KOFFiGIdf","KOFEcGIdj","KOFEcGIdf","IMFnn","spell","SCsEXT","SCsFIN","BA2TOT","cBATOT","BA2FP","SCB","QCB","SCsTOT","QCsTOT","PAsTOT","fi_index_cl_QOG")]
  
  chw=read.dta("E:\\_imf\\CBICND\\Data\\Chwieroth_data.dta")
  chw=chw[,c("iso3","year","areausuk","propteamB","interC")]
  
  kern=read.dta("E:\\_imf\\IMFOWN\\Analysis_rescaled\\IMFCBI_controls.dta")
  kern$lntroops=log(1+kern$ustroops)
  kern$ustroops=NULL
  
  rev=read.dta("E:\\_imf\\IMFOWN\\Analysis\\INETReviewNumbers.dta")
  rev.collapsed=as.data.frame(rev %>% group_by(iso3, year) %>% summarise(events=max(events), maxdur=max(maxdur), maxamou=max(maxamou), reneg=max(reneg)))
  prgtype=read.dta("E:\\_imf\\IMFCMP\\Analysis_Rev\\IMFCMP data_rev.dta")
  prgtype=prgtype[,c("iso3","ayear","a_ebm","a_type")]
  prgfail=read.dta("E:\\_imf\\IMFCMP\\Analysis\\IMFCMP data_20190330.dta")
  prgfail=merge(prgfail,prgtype, by=c("iso3","ayear","a_ebm"), all.x=T, all.y=F)
  prgfail=prgfail[,c("iso3","ayear","inttemp","intpermcor","waiv2","totb","adur","aamou","loanquota0","a_type")]
  prgfail$swaiv=prgfail$waiv2/prgfail$totb*100
  
  
  d=merge(d, chw, by=c("iso3","year"), all.x=T, all.y=F)
  d=merge(d, kern, by=c("iso3","year"), all.x=T, all.y=F)
  d=merge(d, rev.collapsed, by=c("iso3","year"), all.x=T, all.y=F)
  #d=merge(d, prgfail, by.x=c("iso3","year"), by.y=c("iso3","ayear"), all.x=T, all.y=F)
  
  library(haven)
  setwd("C:\\berni\\forschung\\imfown\\submissions\\isq")
  write_dta(d, "IMF spells.dta", version=11)
  
  
  